/*
 * PrimariesGenerator.cc
 *
 *  Created on: 13.02.2011
 *      Author: ed
 */

#include "PrimariesGenerator.hh"
#include "Randomize.hh"

#include "G4Event.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4OpticalPhoton.hh"
#include "G4ParticleMomentum.hh"
G4int PrimariesGenerator::particlesPerEvent;
G4double PrimariesGenerator::photonEnergy = 1*eV;
PrimariesGenerator::PrimariesGenerator(G4int NParticlesPerEvent=500)
{
	particlesPerEvent=NParticlesPerEvent;
	particleGun = new G4ParticleGun(particlesPerEvent);
	particleGun->SetParticleDefinition(G4OpticalPhoton::OpticalPhotonDefinition());
	particleGun->SetParticlePosition(G4ThreeVector(9.0*cm,4.0*cm,9.0*cm));
	SetOptPhotonPolar();

	particleGun->SetParticleEnergy(photonEnergy);

	particleGun->SetParticleMomentumDirection(G4ThreeVector(-1,0,0));

}

PrimariesGenerator::~PrimariesGenerator()
{
	delete particleGun;
}
void PrimariesGenerator::GeneratePrimaries(G4Event* anEvent)
{
	SetOptPhotonPolar();
	particleGun->GeneratePrimaryVertex(anEvent);
}
void PrimariesGenerator::SetOptPhotonPolar()
{
	G4double angle = G4UniformRand() * 360.0*deg;
	SetOptPhotonPolar(angle);
}



void PrimariesGenerator::SetOptPhotonPolar(G4double angle)
{
	if (particleGun->GetParticleDefinition()->GetParticleName() != "opticalphoton")
	{
		G4cout << "--> warning from PrimaryGeneratorAction::SetOptPhotonPolar() :"
			   "the particleGun is not an opticalphoton" << G4endl;
		return;
	}

	G4ThreeVector normal (1., 0., 0.);
	G4ThreeVector kphoton = particleGun->GetParticleMomentumDirection();
	G4ThreeVector product = normal.cross(kphoton);
	G4double modul2       = product*product;

	G4ThreeVector e_perpend (0., 0., 1.);
	if (modul2 > 0.) e_perpend = (1./std::sqrt(modul2))*product;
	G4ThreeVector e_paralle    = e_perpend.cross(kphoton);

	G4ThreeVector polar = std::cos(angle)*e_paralle + std::sin(angle)*e_perpend;
	particleGun->SetParticlePolarization(polar);
}
